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Abstract 

In this paper, an efficient three-dimensional lattice Boltzmann (LB) model with multiple- 
relaxation-time (MRT) collision operator is developed for the simulation of multiphase flows. This 
model is an extension of our previous two-dimensional model (H. Liang, B. C. Shi, Z. L. Guo, and Z. 
H. Chai, Phys. Rev. E. 89, 053320 (2014)) to the three dimensions using the D3Q7 (seven discrete 
velocities in three dimensions) lattice for the Chan-Hilliard equation (CHE) and the D3Q15 lat¬ 
tice for the Navier-Stokes equations (NSEs). Due to the smaller lattice-velocity numbers used, the 
computional efficiency can be significantly improved in simulating real three-dimensional flows, and 
simultaneously the present model can recover to the CHE and NSEs correctly through the chapman- 
Enskog procedure. We compare the present MRT model with the single-relaxation-time model and 
the previous three-dimensional LB model using two benchmark interface-tracking problems, and 
numerical results show that the present MRT model can achieve a significant improvement in the 
accuracy and stability of the interface capturing. The developed model is also able to deal with 
multiphase fluids with very low viscosities due to the using of the MRT collision model, which 
is demonstrated by the simulation of the classical Rayleigh-Taylor instability at various Reynolds 
numbers. The maximum Reynolds number considered in this work reaches up to 4000, which is 
larger than those of almost previous simulations. It is found that the instabilty induces a more 
complex structure of the interface at a high Reynolds number. 
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I. INTRODUCTION 


In the pase two decades, the lattice Boltzmann (LB) method has received great suc¬ 
cess in modeling various fluid systems, and the capability in simulating multiphase flows 
can be recognized as its unique advantage that distinguishes from the traditional numeri¬ 
cal methods [1]. One of most fundamental problems in LB method for multiphase flows is 
how to describe the interfacial dynamics, which is a natural consequence of intermolecular 
interactions between different phases. Up to now, based on different physical pictures of the 
interactions, several types of LB multiphase models has been established, which are substan¬ 
tially divided into four categories: the colour model [2], the pseudo-potential model [3, 4], 
the free-energy model [5, 6] and the phase-field based model [7-11]. In most of these LB 
multiphase models [2-6], the interface is not tracked explicitly, and the region with non-zero 
density gradient is identified as the interface. Therefore the physics of interface tracking 
equation is unknown. Fortunately, the phase-field theory provides a firm foundation on the 
interface physics, in which the interface is tracked by an order parameter that mimics the 
Cahn-Hilliard equation (CHE) [12, 13], 

^ + V-0u = V-M(V/r), (1) 

where cj) represents an order variable, M is the mobility coefficient, /.i is the chemical potential 
and a function of 0, 

/i = 4/50(0 — l)(</> + 1) — kV 2 (p, (2) 

where f3 and k are related to the interface thickness D and surface tension a by the rela¬ 
tionships k = 3 Da/8 and f3 = 3cr/4D . u in Eq.(l) is the fluid velocity and governed by the 
incompressible Navier-Stokes equations [12, 13] 

V • u = 0, (3a) 

p(— + u • Vu) = -Vp + V • [up(V u + Vu T )] + F s + G, (3b) 

where p is the fluid density, p is the pressure, u is the kinematic viscosity, F s is the surface 
tension, and G is the external force. 

Some researchers have constructed some LB multiphase models based on the phase-field 
theory, where the interface is needed to be tracked explicitly by an index or order distribu¬ 
tion function [7, 10, 11, 14], He et al. [7] proposed a LB model for incompressible multiphase 
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flows, in which they adopted an index distribution function to track the interface and a pres¬ 
sure distribution function for solving the flow held. Based on this model, they successfully 
simulated the two-dimensional Rayleigh-Taylor instability, and later the three-dimensional 
case using D3Q15 lattice structure in LB equations for both the interface capturing and 
how held [15]. Although this model is rather robust, it suffers from some limitations, one 
of which is that the recovered interface equation is inconsistent with the CHE noticed by 
Zheng et al. [14, 16]. To this end, they developed a LB model for the CHE, in which a source 
term on a spatial difference of the distribution function is introduced [16]. The model is 
subsequently extended to the three-dimension using D3Q7 lattice model [17]. Recently, Zu 
et al. [10] introduced another similar LB model for the CHE where a spatial difference term 
on the equilibrium distribution function was included. They also modihed the equilibrium 
distribution function in LB equation for flow held such that the continuity equation (3a) 
can be derived. However, the computation of the macroscopic pressure and velocity in their 
scheme is implicit. More recently, we proposed a novel LB model for two-dimensional multi¬ 
phase hows [11]. On one hand, A simpler time-dependent source term is incorporated in LB 
equation for the interface capturing. As a result, the two-dimensional CHE can be recov¬ 
ered correctly. On the other hand, a equilibrium distribution function is delicately designed 
for how held to derive the correct continuity equation while the hydrodynamic properties 
can be computed explicitly. This improved model is also extended to study axisymmetric 
multiphase hows [18] 

As a continuing work, in this paper an efficient three-dimensional LB model for incom¬ 
pressible multiphase how systems is developed based on the multiple-relaxation-time (MRT) 
method. This model has some distinct advantages. Firstly, the model for the CHE requires 
only seven discrete velocities in three dimensions (D3Q7), therefore the expenditure in data 
storage and computational time is smaller than that of other models using D3Q15 lat¬ 
tice [10, 15]. Secondly, the MRT collision model is adopted, which has a better accuracy 
and stability than the single-relaxation-time (SRT) model used commonly in other LB mul¬ 
tiphase models [8-10, 14, 15, 17]. Finally, the present model is able to deal with fluid hows 
at a large Peclet number or a high Reynolds number. The rest of this paper is organized as 
follows. Sec. II presents our three-dimensional MRT LB model for multiphase hows. The 
model then is verihed by several classical numerical experiments in Sec. III. Finally, we 
made a brief summary in Sec. IV. 
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II. THREE-DIMENSIONAL MRT LB MODEL FOR MULTIPHASE FLOWS 


A. Three-dimensional MRT LB model for the Cahn-Hilliard equation 


The LB equation with a MRT collision model for the CHE can be written as, 

/*(x + c iSt, t + 5t) - /i(x, t) = —(M _1 S / M)jj [fj (x, t) - f)] + <5 t Fi(x, t), (4) 

where /j(x, t) is the order distribution function used to track the interface, /® 9 (x, t) is the 
equilibrium distribution function defined as [11, 19] 

0 + (a;* - l)?7/i, i = 0 


/r(x,^) = 


Uji’qu + Ui^- P, i^O, 


( 5 ) 


where the discrete velocities c*, the weighting coefficients a;, and the sound speed c s depend 
on the choice of the discrete-velocity model, 77 is a parameter related to the mobility. In this 
work, an efficient D3Q7 discrete-velocity model is adopted for the CHE, and c., then can be 
given by 

I > 

01-10 0 0 0 

0 0 0 1 -1 0 0 ), ( 6 ) 

0 0 0 0 0 1 -1 

and further according to the following equations (I is the unit matrix) 

^ 1, ^ ' hJjCj 0, (UjCjCj C S I, (7) 

i i i 

u>i can be obtained as 


— 1 — 3A, Wi_6 — —, 


( 8 ) 


where A is a free parameter, c 2 = Ac 2 . To ensure the positive weighting coefficients, the 
parameter A should satisfy 0 < A < The transformation matrix M of the D3Q7 model 
is defined by [ 20 ] 

(\ 1 1 1 1 1 1 \ 


M = 


0 1 -1 0 0 0 0 

0 0 0 1 -1 0 0 

000001-1 
6 -1 -1 -1 -1 -1 -1 

0 2 2 -1 -1 -1 -1 

\0 0 0 1 1 -1 -1 


( 9 ) 
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which is constructed based on the polynomial set of the discrete velocities. in Eq. (1) is 
a diagonal relaxation matrix, 

S f = diag(s 0 ,s 1 , ...,s 6 ), (10) 


where 0 < + < 2, and if the parameters s t equal to each other, the MRT model can reduce 
to the SRT model. To recover the correct CHE, the source term Fj in Eq. (9) should be 
defined as 


F, = [M-'(I - t )M] aRj, 


( 11 ) 


where Ri is given by [11] 


Ri = 


UiCi ■ d t (j )u 


In the present model, the order parameter is computed by 


( 12 ) 


4> = F,f" 

i 

and the density is taken as a linear function of the order parameter, 

_ 1 + 0 , 1-0 
P ~ 9 Pi - Pgi 


(13) 


(14) 


where pi and p g represent the densities of liquid and gas phases, respectively. 

The evolution of LB equation (4) can be commonly divided into two steps, i.e., the 
collision process, 


ft = /i(x,t) - (M^S'M)*- [/j(x, t) - /f (x,*)] + 5 t Fi(x,t), (15) 

and the propagation process, 


/j(x + Cj St,t + St) = . 


(16) 


To reduce the matrix operations, it is wise that the collision process of MRT model is 
implemented in the moment space. By premultiplying the transformation matrix, we can 
easily derive the equilibrium distribution function in moment space, 


mf"' = (0, 60 _ 2L4///J, 0,0) T , 


(17) 


c c c 

where u x , u y and u z are the x-, y- and z- components of macroscopic velocity u, respectively. 
Similarly, the source term Ri in the moment space can be presented as 


mR = (o, 0, 0, 0) T 


(18) 
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The Chapman-Enskog analysis is carried out on the LB evolution equation (4), and the 
results demonstrate that the CHE can be derived correctly from the present MRT model, 
and the relationship between the mobility M and the relaxation parameter is also derived 
as 

M = r]c 2 s (Tf — 0.5)ht, (19) 

where t/ = 1 /si, and Si = S 2 = S 3 . 


B. Three-dimensional MRT LB model for the Navier-Stokes equations 

The MRT LB equation with a source term for the NSEs reads as [21] 


gi(x + ciS t , t + S t ) - gi(x, t ) 


-(T 1 S s r) ij [ 5 f i (x,t) - gf(x,t)] + StGi, 


( 20 ) 


where g t is the density distribution function, g\ q is the equilibrium distribution function and 
is defined by 


9i 


eq — 


fs(u>i - 1) + pSi(u) i = 0 
%Ui + p8i( u) i^O 


( 21 ) 


with 


Sj(u) = LOi 



+ 


(Cj •u ) 2 

2 ci 


u • u 

"2cT 


( 22 ) 


To simulate the fluid flows in three dimension, we can choice several types of lattice velocity 
models, such as D3Q15 or D3Q19 [ 22 ]. The D3Q15 lattice model is used in this work due 
to its smaller data storage and higher computational efficiency. Following the work of Qian 
et al. [22], the discrete velocities c., of the D3Q15 lattice model can be given by 


/ 


c i = c< 


V 


0 1 -1 0 0 0 0 1 
00 0 1-10 0 1 
000001-11 


-1 

-1 

-1 


\ 

1-11-11-1 
1 —1—1 1 —1 1 >, 
-11 1-1-11 


(23) 


where c = \/3 c s . According to the ordering of c,, the weight coefficients are presented as 


_ 2 _ 1 1 
M) — g, kfi_6 — -, W7-14 — 72 ' 


( 24 ) 
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The transformation matrix T in Eq. (20) can be given by [23] 
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-1 
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-1 

-1 

1 

1 

1 
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1 

1 
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1 

16 
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-4 

-4 
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-1 
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1 
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and the corresponding diagonal relaxation matrix S 9 is denoted by 


(25) 


S 9 = diag(X 0 ,X 1 , ...,X 14 ), 


( 26 ) 


where 0 < A* < 2. The source term Gi in Eq. is defined as [11] 

Qfl 

Gi = [T-‘(I - y)r]yTj 


( 27 ) 


where 


Ti — -- 2 -“ ' [ s i( u )V(p c s) + (F s + F a + G)(Sj(u) + Wj)] , 


( 28 ) 


in which F a = 0.5(p4 — pb)M'V 2 pu is an interfacial force, and F s is the surface tension 
taken the potential form F s = p’V(f). In the present model, the macroscopic pressure and 
velocity can be obtained from 

X! c i9i + 0.55f(F s + G) 

( 29 ) 


u = 


p - 0.25 (p A - pb)M\7 2 p 


V 


(1 — Wo) 


+ Vp + ps°(u) 

Li^o 


( 30 ) 
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The collision process of MRT LB equation for NSEs is also carried out in the moment 
space. After some algebraic manipulations, the equilibrium distribution function g f - q and 

source term T) in the moment space can be respectively derived as 

3p + pn 2 45p + 5pu 2 pu x 7pu x pu y 7pu y pu z 7pu z 


mg'" = (0, 


O) 5 5 o ? 5 o ’ 5 O ? 

c z c c Sc c Sc c Sc 

^ u x ^z ^y ^ z P^x^y pUyll z pU x li z fq-\\ 

P o ‘iP o’ o > “2 1 “2 5 


and 


mT = [u • Vp, 


u-(-Vpc 2 + 2F) u • (7Vpc 2 + 10F) F x 7F X F y 7F y F z 7F z 

c z c z c Sc c Sc c Sc 

Au x (d x pc 2 s + F x ) - 2 u y (d y pc 2 s + F y ) - 2 u z (d z pc 2 s + F z ) 2 u y (d y pc 2 s + F y ) - 2 u z (d z pc 2 s + F z ) 

c 2 ’ c 2 

U x (d y pc 2 + Fy) + u y (d x pc 2 + F x ) u y (d z pc 2 + F z ) + u z (d y pc 2 + F y ) 


HxidzPCg T Fz) T l U'z(d x pc 8 T -Fx) „n r p 

12 ’ 

c z 

(32) 

We also conduct the Chapman-Enskog analysis on LB Equation (20), and the results 
show that the present model can correctly recover to the NSEs with the kinematic viscosity 
determined by 

v = c s( T <? “ 0.5)5 t , (33) 

where t 9 = 1/A 9 , and A 9 = Ai 0 = An = X u = A i3 . 

In practical applications, the derivative terms in present model should be discretized 
by suitable difference schemes. As widely adopted in the references [19, 24], an explicit 
difference scheme 

x(x,t) - xCM- St) 


<%x(x, t) = 


(34) 


is used for computing the time derivative in Eq. (12), and the isotropic central schemes [25] 

w,qx(xbq(5 t ,t) 


Vx( x , t) = 

i^O 


c 2 Jt 


and 


V 2 x(x,t) = Y, 

0 


2uiiCi[x(x+ c,S,,t) - x(x,«)l 




(35) 


(36) 


are employed for calculating the gradient and the Laplacian operator, respectively. In above 
equations, x represents an arbitrary function. It should be noted that the schemes (35) 
and (36) not only can preserve a secondary-order accuracy in space, but also can ensure the 
global mass conservation of a multiphase system [25]. 
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(b) 

FIG. 1. The rotation of Zalesak’s sphere during one period at Pe = 200: (a) the present MRT 
model; (b) the previous three-dimensional LB model [17]. 

III. NUMERICAL RESULTS AND DISCUSSIONS 

In this section, we will validate the present three-dimensional MRT LB model with several 
numerical examples. We first test the performance of the three-dimensional LB model in 
interface capturing by simulating two classical benchmark problems: rotation of the Zalesak’s 
sphere [26] and deformation field flow [27]. In these simulations, the evolution equation (4) 
is only adopted in that the velocity distribution has been specified in advance. Next, we 
simulated the three-dimensional Rayleigh-Taylor instability to show the ability of the present 
model for multiphase flows, where a comparison between the numerical results and some 
available results is also conducted. 



A. Rotation of the Zalesak’s sphere 

Rotation has been widely used in the literatures to test interface tracking methods [11, 
26, 28]. Here we consider the rotation of the Zalesak’s sphere which has a slot of 16 lattice 
units width in a 100 x 100 x 100 domain. This sphere is initially centered at (50, 50, 50) and 
the radius R occupies 40 lattice units. The revolution of the sphere is driven by a constant 
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FIG. 2. The rotation of Zalesak’s sphere during one period at Pe = 600: (a) the present MRT 
model; (b) the present SRT model. 


vorticity velocity field, 


u(x,y,z) = 0, 

v(x, y, z) = ([/7t/100)(50 — z), (37) 

w(x, y, z) = (C/7r/100)(y - 50), 

where U takes a value of 0.02 so that the sphere complete one cycle every 10 4 time steps. In 
our simulations, the interface thickness D and surface tension a are fixed as 2.0 and 0.04, 
respectively; the relaxation matrix is set as 

S / = d^(1.0 ,^,|l.2 ,1.2,1.2). (38) 


The periodic boundary condition is applied at all boundaries. To examine the mobility 
effect, we introduce the dimensionless Peclet number, which is defined as, 


Pe = 


LU 

H' 


(39) 


where L is the characteristic length, and set as the value of R. Different mobilities can be 
obtained by changing the value of Pe. Figure 1 shows the revolution processes of the sphere 
during one period at Pe = 200. It is seen from Fig. 1 that the present model can precisely 
capture the interface shape in one period, and the sphere returns to its initial configuration 
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FIG. 3. Comparison of the relative errors produced by different LB models in Zalesak’s sphere 
rotation tests with various Peclet numbers. 

at time T, which is accordance with the expected results. For a comparison, the results 
obtained by the previous three-dimensional LB model [17] are also presented in Fig. 1. It is 
observed that the previous model not only produces some obvious sawteeth on the sphere 
surface, but also induces some unphysical disturbances around the computational domain. 
This implies that the present model can obtain a more accurate and stable interface. It has 
been mentioned above that the SRT model can be deemed as a special case of the MRT 
model when all the relaxation parameters equal to each other. Due to the more adjustable 
relaxation factors, the MRT model should show more potential to achieve a better numerical 
stability against the SRT model. To illustrate this point, we present a comparison between 
them using this case. Figure 2 depicts snapshots of rotating sphere during one period at 
Pe = 600 by the MRT and SRT models. It is clearly seen that the results of the SRT model 
are unstable. The slot of the sphere is slightly distorted and there produces some extra 
jetsam at the corners of the computational domain. In contrast, the present MRT model 
can capture the moving interface of the sphere correctly. To further quantitatively describe 
the accuracy of the interface capturing, the global relative error on the order parameter is 
introduced as 

Frr = Ex|0(x,r)-0(x,o)| 

Ex|0(x,o)| 



li 
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\t 





T 





y| t 



\l \t 
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(c) 


FIG. 4. Snapshots of deformation field test at Pe = 400: (a) the present MRT model; (b) the 
present SRT model; (c) the previous LB model [17]. The times from the left pattern to the right 
are 0, T/4, T/2, 3T/4 and T. 


where </>(x, T) is the value of the order parameter at position x and time T, and </>(x, 0) is 
the exact solution at initial time. We computed the relative errors generated by the LB 
previous model [17] and the present MRT and BGK models at various Peclet numbers, and 
show the results in Fig. 3. As seen from this figure, one can find that the MRT model is 
more accurate than the SRT model, which is further more accurate than the previous model, 
especially at large Pelcet number. The effect of the model parameter A is also examined, 
and it is found from Fig. 3 that the parameter A has almost no effect on the accuracy of the 
MRT model. Without loss of generality, in the following simulations we fix the parameter 
A at 0.25. 
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-♦-Present MRT model 
—b— Present SRT model 



0 50 100 150 200 250 300 350 400 450 500 

Pe 


FIG. 5. Comparison of the relative errors produced by different LBMs in deformation field tests 
with various Peclet numbers. 

B. Deformation field flow 

The above test does not induce large changes of the interface. To show the capacity of 
the present MRT model, in this subsection we will consider a rather challenging problem 
of deformation held flow, in which the interface could undergo a large deformation [26, 27]. 
The initial setup of this problem is described as follows. A sphere with a radius of 15 is 
placed in a 100 x 100 x 100 computational domain centered at (35, 35, 35). The velocity 
held of this flow is strongly nonlinear and is given by 

u(x,y,z ) = 2Usin 2 (7nr/100)sin(27n//T00)sin(27r,2:/100), 

v(x,y,z ) = — Usin(27nr/100)sin 2 (7n//100)sin(27r,2/T00), (41) 

w(x,y,z) = — Tsin(27r:r/100)sin(27n//100)sin 2 (7r;2/100), 

and we postmultiply a time-dependent function cos(7 rt/T) to make the flow periodic, where 
t is the iteration step divided by 100/17. In our simulations, U and T are given as 0.02 and 
2; the other physical parameters and boundary condition are set as those in the last test. In 
theory, the sphere will undergo the deformation continuously until time T/2 , and the velocity 
field then is revised in time, the sphere goes back and returns to its original position. Figure 
4 shows the evolution of the interface pattern at Pe = 400 obtained by the use of the MRT 
model, the SRT model and the previous LB model [17]. It can be clearly observed that the 
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present MRT model can give an accurate prediction in the evolution of the interface: it has 
a largest deformation at time Tj 2 and moves back to the initial configuration at time T. 
The behaviors of the interface conform to the theoretical results. In contrast, the results of 
the SRT model and the previous three-dimensional LB model are unstable. At time T/4 
or T, some jagged shapes are produced by the SRT model in the vicinity of the interface. 
At the same time, the previous LB model is worst in tracking the interface. A massive 
amount of unphysical disturbances can be clearly observed in the system. We also perform 
simulations with different Peclet numbers, and compare the relative errors generated by 
these LB models. The results are presented in Fig. 5. From this figure, one can find that 
the present MRT model is more accurate than the SRT model and the previous LB model. 


C. Rayleigh-Taylor instability 

At last, we simulated a benchmark problem of the Rayleigh-Taylor instability (RTI). 
RTI is a classical and common instability phenomenon, which occurs at perturbed interface 
between two different fluids, where the gravity force is applied. RTI studies are useful since 
it has particular relevance and importance in fields including inertial confinement fusion [29] 
and astrophysics [30]. For this reason, RTI has become a subject of intensive researches 
in the past 60 years using theoretical analysis [31], experimental methods [32, 33] and also 
numerical approaches [7,11, 13, 15, 34-37]. However, to our best knowledge, most of previous 
numerical studies are limited to the two-dimensional case [7, 11, 13, 36, 37], and relatively 
few attention has been paid on the three-dimensional RTI [15, 34, 35]. On the other hand, 
the simulated Reynolds numbers considered in previous studies are small in general. In this 
subsection, we applied the present MRT model to study the three-dimensional Rayleigh- 
Taylor instability, and the effect of Reynolds number on the evolution of the interface was 
examined. 

The physical problem we considered here is a rectangular box with an aspect ratio of 
4A x A x A, where A is the box width. The initial interface is located at the midplane 
(z = 2A), with an imposed square-mode perturbation 

h(x, y) = 0.05A[cos(—+ cos(—y^)], (42) 

A A 
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(a) 


(b) 


FIG. 6. Evolution of the fluid interface in immiscible RTI at different Reynolds numbers: (a) 
Re=1024; (b) Re=4000. The corresponding times from the left pattern to the right are t=1.0, 2.0, 
3.0, 4.0, where t is normalized by the characteristic time X/^/~g\. 


and the initial order distribution then can be given by 

2 [z-h- 2A) 


<j){x ) y , z) — tanh ■ 


D 


According to [15], the Reynolds numbers (Re) characterizing RTI can be defined as 

\fg\\ 


Re — 


v 


(43) 


(44) 


where g is a gravitational acceleration, and v is the kinematic viscosity. In our simulations, 
we take the densities of liquid and gas phases as 3.0 and 1.0, corresponding to an Atwood 
number of 0.5; some physical parameters are given as: A = 128, \fg\ = 0.04, A = 128 
D = 4.0, a = 0.0001; the Peclet number defined in [10,11] is fixed at 50.0; the relaxation 
matrix is set to be 

S / = diag( 1.0,1.25,1.25,1.25,1.2,1.0,1.0), (45) 


and the relaxation factors in S 9 are chosen to be unity expect for the ones related to the 
kinematic viscosity. The periodic boundary conditions in the lateral directions and no-slip 
boundary condition in the vertical direction are applied in our studies. Figure 6(a) shows the 
evolution of the density contours in immiscible RTI at a Reynolds number of 1024. It can be 
seen that, due to the gravity effect, a heavy fluid and a light one penetrate into each other at 
early time, and then forms the spike and bubble, respectively. After that, the heavy fluid rolls 
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up along the flank of the spike, and a mushroom-like structure appears (t=3.0), which can 
be attributed to Kelvin-Helmholtz instabilities providing the rolling motion of the interface. 
Finally, the mushroom develops further and becomes much bigger. The patterns of the liquid 
interface obtained by the present model compare well with previous results [10, 15, 34], We 
also simulated the three-dimensional RTI at a high Reynolds number of 4000, which has not 
been considered in Ref. [15], and presented the results in Fig. 6(b). It can be observed that 
the interface takes on almost similar behaviors before time 3.0. However, it subsequently 
presents a distinct manner. Due to larger shear interaction between different layers, the 
interface becomes more complex, which eventually undergoes a breakup inducing some tiny 
dissociative drops in the system. To observe the evolution of the interface more clearly, we 
plotted in Fig. 7 the interface patterns at the diagonal vertical plane (x = y ) with above two 
Reynolds numbers. It is found that for both Re = 1024 and Re = 4000, the development 
of the initial mode follows the pattern similar from two-dimensional simulations [11]. In the 
following, two pairs of counter-rotating vortices are formed at the spike tip and the saddle 
point (see time 3.0), which is significantly different from the two-dimensional results [11]. 
The vortices grow with time for Re = 1024, while they become unstable for Re = 4000, 
resulting in the mixing of two different fluids at the vicinity of the spike tip. We also gave 
a quantitative study on the Reynolds number effect. Figure 8 depicts the evolution of the 
positions of the bubble front and spike tip obtained by the present MRT model and the 
previous numerical results in Ref. [15]. It is shown that the results of Re = 1024 obtained 
by our model agree well with those in Ref. [15], which verifies the numerical accuracy of the 
present MRT model in dealing with complex interfacial flows. From Fig. 8, one can also find 
that there is no evident difference in trajectory of bubble front at two different Reynolds 
numbers, while the spike tip moves slightly faster at a larger Reynolds number. 


IV. SUMMARY 

In this work, an efficient three-dimensional lattice Boltzmann model based on the MRT 
collision method is proposed for multiphase flow systems. This model is a straightforward 
extension of our previous two-dimensional model [11] to the three-dimensions. The present 
model for the CHE only utilizes seven discrete velocities in three dimensions, while most of 
previous LB models need at least fifteen discrete velocities. As a result, the computational 
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FIG. 7. Evolution of the fluid interface at the diagonal vertical plane (x — y ) with different 
Reynolds numbers: (a) Re=1024; (b) Re=4000. The corresponding times from the left pattern to 
the right are t=1.0, 2.0, 3.0, 4.0, where t is normalized by the characteristic time A /yfg\. 



FIG. 8. The time evolution of the positions of the bubble front, spike tip. The length and time 
are normalized by A and A/y^A, respectively. 

efficiency of the present model can be greatly improved in simulating three-dimensional mul¬ 
tiphase flows. In addition, the advanced MRT collision model is adopted in LB equations for 
both the CHE and the NSEs, which has a better stability than the SRT model. Two classical 
interface-capturing problems including rotation of the Zalesak’s sphere and deformation field 
flow were conducted to test the model, and the results show that the present MRT model 


17 











is more stable and accurate than the SRT model and the previous LB model in tracking 
the interface. Finally, the present MRT model is applied to study the three-dimensional 
Rayleigh-Taylor instability at various Reynolds numbers. It is found that the numerical 
results at low Reynolds numbers agree well with previous data, while the instability at a 
high Reynolds numbers induces a more complex structure of the interface. 
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